MODULE ZVAR
!
! Verwendungszweck:
!
!      Bereitstellung von Hilfsvariablen fuer die Integration der
!      hydrostatischen Schichtung in z-Richtung
!
!-----------------------------------------------------------------------
! Interne Variable:
!
!     pseudo2D : .true.: z-integrierte Dichte und Hydro-gleichungen, .false.: nur in z=0 Ebene
!
!     H_p      : Druckskalenhohe der hydrostatischen Schichtung in z-Richtung, neuer Zeitpunkt
!     H_pZ     : Druckskalenhohe der hydrostatischen Schichtung in z-Richtung, zeitlich zentriert
!     dH_pZdE  : Ableitung von H_pZ nach der inneren Energie (zum neuen Zeitpunkt)
!     dH_pZdR  : Ableitung von H_pZ nach dem Radius (zum neuen Zeitpunkt)
!
!     rho0     : Dichte in der Mittelebene der Scheibe, d.h. fuer z=0
!
!     H_p_flag : Schalter fuer die Neigung der Schabenoberflaeche, fuer Einstrahlung in Eq_ene
!
!     z_gamma  : Adiabatenexponent, fuer die Integration der hydrostatischen Schichtung in z-Richtung
!                1-atomige Gase:                            c_v = 3/2, c_p = 5/2, gamma = 5/3
!                2-atomige Gase:                            c_v = 5/2, c_p = 7/2, gamma = 7/5
!                3-atomige Gase, starre Atombindung:        c_v = 6/2, c_p = 8/2, gamma = 8/6
!                3-atomige Gase, nicht starre Atombindung:  c_v = 7/2, c_p = 9/2, gamma = 9/7
!
!-----------------------------------------------------------------------

      use config, only : rkind, np
      use physco, only : z0

      implicit none

      private
      public :: ZINT, pseudo2D, H_p, dH_pdE, dH_pdR, H_pZ, dH_pZdE, dH_pZdR, rho0, rho0Z, H_pZ_flag, z_gamma


      logical          :: pseudo2D

      real(kind=rkind) :: H_p(np),  dH_pdE(np),  dH_pdR(np)
      real(kind=rkind) :: H_pZ(np), dH_pZdE(np), dH_pZdR(np)
      real(kind=rkind) :: rho0(np), rho0Z(np)
      real(kind=rkind) :: H_pZ_flag(np) = z0

      real(kind=rkind) :: z_gamma

CONTAINS





SUBROUTINE ZINT
!
! Verwendungszweck:
!
!     Berechnung der Druckskalenhohe der hydrostatischen Schichtung in z-Richtung
!     und Umrechnung der z-integrierten Dichte in die Dichte fuer z=0
!

      use primvar, only : X, XZ, MR, MD, ME
      use config,  only : np
      use global,  only : zz
      use physco,  only : grav, sqrt2pi, z0, z1, z3, z12
      use RBvar,   only : M_star, H_star
      use RBvar,   only : R_in, rho0_in, rho_in, e_in, R_out, rho0_out, rho_out, e_out


      implicit none

      integer       :: idx
      logical, save :: init = .true.


!-----------------------------------------------------------------------
!    Druckskalenhoehe und Dichte fuer z=0
!-----------------------------------------------------------------------

      if (pseudo2D) then

         H_p(:)     = sqrt(  X(ME,:) * (z_gamma-z1) *  X(MR,:)**3 / (grav * M_star) )
         H_pZ(:)    = sqrt( XZ(ME,:) * (z_gamma-z1) * XZ(MR,:)**3 / (grav * M_star) )


         dH_pdE(:)  = z12 / H_p(:)   *  (z_gamma-z1) *  X(MR,:)**3 / (grav * M_star)
         dH_pZdE(:) = z12 / H_pZ(:)  *  (z_gamma-z1) * XZ(MR,:)**3 / (grav * M_star) * zz

         dH_pdR(:)  = z12 / H_p(:)   *   X(ME,:) * (z_gamma-z1) * z3 *  X(MR,:)**2 / (grav * M_star)
         dH_pZdR(:) = z12 / H_pZ(:)  *  XZ(ME,:) * (z_gamma-z1) * z3 * XZ(MR,:)**2 / (grav * M_star) * zz

      else

         H_p(:)     = z1 / sqrt2pi
         H_pZ(:)    = z1 / sqrt2pi


         dH_pdE(:)  = z0
         dH_pZdE(:) = z0

         dH_pdR(:)  = z0
         dH_pZdR(:) = z0

      end if

      rho0(:)  = X(MD,:)  / ( sqrt2pi * H_p(:)  )
      rho0Z(:) = XZ(MD,:) / ( sqrt2pi * H_pZ(:) )


! Schalter fuer die Einstrahlung
      do idx = 3,np-3   ! hier kann man den physikalich berenzt sinnvollen Einstrahlungsterm fuer die aeusserste skalare Zelle, idx = np-3, ausschalten
         if ( (H_pZ(idx+1) - H_star) / XZ(MR,idx+1) > (H_pZ(idx) - H_star) / XZ(MR,idx) ) then
            H_pZ_flag(idx) = z1
         else
            H_pZ_flag(idx) = z0
         end if
      end do



!-----------------------------------------------------------------------
!    Integrierte Saeulendichte fuer die Randbedingungen
!-----------------------------------------------------------------------

      if (init) then
         if (pseudo2D) then

            rho_in  = rho0_in  * sqrt2pi * sqrt( e_in  * (z_gamma-z1) * R_in**3  / (grav * M_star) )
            rho_out = rho0_out * sqrt2pi * sqrt( e_out * (z_gamma-z1) * R_out**3 / (grav * M_star) )

         else

            rho_in  = rho0_in
            rho_out = rho0_out

         end if
         init = .false.
      end if


END SUBROUTINE ZINT


END MODULE ZVAR




